Anisotropic total variation minimization approach in in-line phase-contrast tomography and its application to correction of ring artifacts
Ji Dong-Jiang1, 2, Qu Gang-Rong1, 3, †, Hu Chun-Hong4, Liu Bao-Dong5, 6, 7, Jian Jian-Bo8, Guo Xiao-Kun9
School of Science, Beijing Jiaotong University, Beijing 100044, China
School of Science, Tianjin University of Technology and Education, Tianjin 300222, China
Beijing Center for Mathematics and Information Interdisciplinary Sciences, Beijing 100048, China
College of Biomedical Engineering, Tianjin Medical University, Tianjin 300070, China
Division of Nuclear Technology and Applications, Institute of High Energy Physics (IHEP), Chinese Academy of Sciences (CAS), Beijing 100049, China
Beijing Engineering Research Center of Radiographic Techniques and Equipment, Beijing 100049, China
University of Chinese Academy of Sciences, Beijing 100049, China
Radiation Oncology Department, Tianjin Medical University General Hospital, Tianjin 300052, China
The Second Hospital of Tianjin Medical University, Tianjin 300070, China

 

† Corresponding author. E-mail: grqu@bjtu.edu.cn

Project supported by the National Natural Science Foundation of China (Grant Nos. 61671004, 61271012, 81371549, 81671683, and 11501415), the Natural Science Foundation of Tianjin City, China (Grant No. 16JCYBJC28600), the WBE Liver Fibrosis Foundation of China (Grant No. CFHPC20131033), the Instrument Developing Project of the Chinese Academy of Sciences (Grant No. YZ201410), the Foundation of Tianjin University of Technology and Education (Grant Nos. KJ11-22 and J10011060321), SRF for ROCS, SEM, and the IHEP-CAS Scientific Research Foundation (Grant No. 2013IHEPYJRC801).

Abstract

In-line phase-contrast computed tomography (IL-PC-CT) imaging is a new physical and biochemical imaging method. IL-PC-CT has advantages compared to absorption CT when imaging soft tissues. In practical applications, ring artifacts which will reduce the image quality are commonly encountered in IL-PC-CT, and numerous correction methods exist to either pre-process the sinogram or post-process the reconstructed image. In this study, we develop an IL-PC-CT reconstruction method based on anisotropic total variation (TV) minimization. Using this method, the ring artifacts are corrected during the reconstruction process. This method is compared with two methods: a sinogram preprocessing correction technique based on wavelet-FFT filter and a reconstruction method based on isotropic TV. The correction results show that the proposed method can reduce visible ring artifacts while preserving the liver section details for real liver section synchrotron data.

1. Introduction

In medical applications, computed tomography (CT) has been widely used to reconstruct the detailed structure of objects. A CT imaging system can be modeled as:

where A denotes the system weight matrix, p denotes the projection data collected by the detector, and v denotes the linear attenuation coefficients of the reconstructed object. The task of CT imaging is to reconstruct v from p, x-ray computed tomography (CT) uses attenuation coefficients of an object to reconstruct interior structures. However, the attenuation values change little between different biological soft tissues, so tissue characteristics are indiscernible in x-ray CT images. In-line phase-contrast (IL-PC) imaging is a new imaging method that generates high spatial resolution and super contrast of weakly absorbing objects compared with conventional attenuation-contrast imaging.[1] In practical biomedical applications, IL-PC is combined with CT, and IL-PC-CT can reconstruct images of soft samples. IL-PC-CT images, reconstructed from the two-dimensional (2D) x-ray projection data acquired with a flat-panel detector, are often corrupted by serious ring artifacts, because of malfunctioning and miscalibration of elements in the x-ray detectors as well as impurities in the scintillator crystals.[2] As visualization and quantification of the anatomic and pathological features are influenced by these ring artifacts, removal or reduction of these artifacts is essential. Ring artifact removal or reduction without impairing the image quality is still a challenging problem for the researchers. To overcome the challenge, several methods have been reported so far to correct the ring artifacts from CT images. The flat field correction is a common method to reduce ring artifacts.[3] However, ring artifacts are not removed completely by this method when different camera elements have intensity dependency, nonlinear response functions or the incident beam has time-dependent non-uniformities.[2] Another method for ring artifacts suppression is based on moving the sample or the detector system during the acquisition in defined horizontal and vertical steps.[4,5] However, this requires very precise translation motor components. The third method for correcting ring artifacts is the image-processing based on sinogram preprocessing and reconstructed images postprocessing. Some reported methods dealing with the raw sinogram aim at detecting and correcting the spurious lines in the sinogram before applying the reconstruction process.[6,7] The post-processing approach may be tempting to detect and remove the ring artifact in the image domain rather than the sinogram domain.[8,9]

TV regularization has been a very active research topic in CT reconstruction.[1013] One example in the field of image reconstruction is the isotropic TV norm.[14,15] Later, the isotropic TV minimization approach was addressed in in-line phase-contrast x-ray tomography.[16] Although the isotropic TV norm was proven effective for low-dose CT and interior tomography, the isotropic TV minimization is unfit for limited-angle CT. A weighted anisotropic TV minimization method[17] is introduced to address the limited-angle CT problem using a simultaneous algebraic reconstruction technique (SART) algorithm[18] and the gradient descent method.[19]

The TV model falls into two categories: isotropic TV and anisotropic TV.[20,21] The mathematical definitions for both isotropic and anisotropic TV are given in the discrete setting. v is denoted as the column vector by a lexicographical ordering of a two-dimensional (2D) image. In our study, the isotropic and anisotropic TV are respectively defined as[22]

where is the first-order finite difference operator in the horizontal and vertical directions.

The TV term (isotropic or anisotropic) is essentially the -norm of the gradient of the image.[23] Thus, the sparsity constraint -norm forces the ring artifacts variables to have only a few not null components, since the -norm is a good approximation of the sparsity-inducing -norm. Similarly, in this study, the ring artifacts can be wide and faint, which can be corrected by the TV term. The isotropic TV equals to the square root of the first-order finite difference operator in the horizontal and vertical directions, which is similar to the -norm of a complex vector.[24] Thus, first-order finite difference operator in the horizontal and vertical directions are unseparated, which leads to the mutual interference. When the images become more complex, i.e., images that contain complex textures and gradual intensity like liver images, isotropic total variation can often produce cartoon-like approximations; for the result see Subsection 3.1. The main idea of anisotropic TV is to separate the first-order finite difference operator in the horizontal and vertical directions to avoid the interference between them, hence the anisotropic total variation regularization can preserve images sharp with corners,[25,26] and natural images usually contain a lot of edges and other anisotropic features, since the anisotropic TV is a good sparsity approximation of natural images. In order to reduce visible ring artifacts while preserve the image section details(edges, corners, and other anisotropic features), based on the advantages of anisotropic TV outlined above, we propose to correct the ring artifacts by use of an anisotropic TV term during the image reconstruction process; for the result see Subsection 3.2.

There are two differences between our method and the weighted anisotropic TV minimization method.[17] Firstly, the weighted anisotropic TV minimization method is introduced to address the limited-angle CT problem, while our aim is to correct image ring artifacts and reconstruct an image from the complete-view without introducing weighting factors for anisotropic TV minimization. Secondly, unlike the calculation method of the weighted anisotropic TV minimization,[17] in our study, the anisotropic TV minimization method is introduced to correct image ring artifacts using a simultaneous algebraic reconstruction technique (SART) algorithm[18] and the alternating direction method of multipliers (ADMM).[27]

The main features and contributions of this work are summarized as follows.

The rest of this paper is organized as follows. In Section 2, we will describe an anisotropic TV-L2 model and its implementation based on the ADMM.[27] In Section 3, the real liver section synchrotron data results are presented for the anisotropic TV-L2 model and demonstrate the effectiveness of the proposed method for correction of ring artifacts. Finally, concluding remarks are provided in Section 4.

2. Method
2.1. In-line phase-contrast imaging

The experiments were performed in BL13W1 of the Shanghai synchrotron radiation facility (SSRF). A liver sample was placed on the shelf, and its image was reconstructed using an IL-PC x-ray imaging setup. Figure 1 shows the schematic of the IL-PC x-ray imaging setup.[28] The monochromator crystal was a combination of an Si (111) orientation crystal and an Si (311) orientation crystal. A high-precision sample platform was used to position the sample and rotate the sample axis perpendicular to the beam. The x-ray beam energy in the experiments was set at 20 keV, and the detector employed an x-ray CCD camera system with per pixel. The projections were recorded with a sample-to-detector distance of 1 m and exposure time of 10 ms (as shown in Fig. 2(a)). In addition, 20 flat field images and 10 dark field images were collected.

Fig. 1. (color online) Schematic diagram of an in-line phase contrast imaging setup.
Fig. 2. (a) Sinogram of the projection data from the liver section:projections used in this sinogram were collected from 180° angle of view (1024 projections with 0.176° projection interval), (b) magnified region of interest (ROI) in the sinogram.

Figure 2(a) shows a sinogram image of the projection data from the liver section. The data were collected by taking 1024 projections with total 180° angle of view. In Fig. 2(b), the magnified ROI of the sinogram reveals vertical lines that are responsible for ring artifacts. In this study, our aim is to correct ring artifacts during the reconstruction process.

2.2. Algebraic reconstruction method

To solve the problem in Eq. (1), we use SART to enforce consistency with the projection data. The formula of SART can be expressed as[18]

where K indicates the iteration number, means the j-th 2D pixel, is the n-th 2D pixel contribution along the i-th ray, is the measured projection, and is the relaxation parameter.

2.3. Anisotropic TV-L2 model

In our study, an image is reconstructed as a solution for the following anisotropic TV-L2 model:

where is the regularization parameter.

A new reconstruction algorithm is developed based on proximal forward–backward splitting (PFBS),[29] which completely decouples data fidelity minimization and image regularization problems. As a result, a CT image reconstruction algorithm step without any regularization operators can be derived for data fidelity minimization, and then sole image regularization or denoising subproblems without projection operators can be rapidly solved using ADMM.[27] The detail of the PFBS algorithm regarding how to decouple data fidelity minimization and image regularization problems was introduced by Guo in 2016.[30] According to the PFBS optimization algorithm, the anisotropic TV-L2 image reconstruction model mainly consists of the following two iterative steps with the initial image :

In Eq. (6), v is the image to be reconstructed, p is the projection data, A the system matrix or a projection operator discretized from the x-ray transform, the step for computing is the gradient descent update with a step size of for the problem . There are many methods to solve the linear inverse problem for v iteratively. In this study, the SART algorithm is used to solve the subproblem Eq. (6).

For the sake of convenience, we can rewrite the subproblem Eq. (7) as

where and .

Next, we provide a solution algorithm for the subproblem Eq. (8) with anisotropic TV regularization Eq. (3) based on ADMM. We adopt a splitting method from Ref. [31] by introducing two auxiliary variables y = f, z = Df, and the optimization is rewritten as

Its corresponding augmented Lagrangian function is

where λ is the regularization parameter, vectors μ and ρ are the Lagrange multipliers, and δ and θ are the penalty parameters. We use the alternating direction method to iteratively solve the following subproblems

The Lagrange multipliers are expressed as

The subproblem is to minimize the augmented Lagrangian function; these subproblems are investigated one by one.

In summary, the implementation steps of the TV-L2 model based on SART-anisotropic TV are given as follows:

2.4. Comparison method
2.4.1. SART based on isotropic TV (SART-isotropic TV)

The CT reconstruction problem given by Eq. (1) can be solved by the compressive sensing-based reconstruction method to minimize the image isotropic TV regularized by the projections. This process is equivalent to solving the following optimization program:

The model described by Eq. (22) was investigated using the alternating minimization procedure from:[14]

The gradient descent is controlled via the step size parameter β. The SART step loop is executed K times. denotes the total number of gradient descents.

2.4.2. SART based on sinogram preprocessing with the wavelet-FFT filter technique(SART-wavelet-FFT filter)

In order to correct ring artifacts, a fast, powerful and stable filter based on combined wavelet and FFT transforms was presented to preprocess the sinogram.[7] The steps of the filter based on combined wavelet and FFT transforms are as follows:

Step 1 in the first step, a tight condensation of the stripe information by combining the wavelet and FFT transforms.

Step 2 it then performs a damping of the relevant coefficients by using a Gaussian function.

When the sinogram preprocessing is complete, the output sinogram can be used to reconstruct the image by the SART algorithm. The reconstruction algorithm based on sinogram preprocessing is referenced as the SART-wavelet-FFT filter.

3. Numerical results

The TV-L2 model based on SART-anisotropic TV can be used to correct the ring artifacts for IL-PC-CT. We compare the algorithm with closely related methods: SART, SART-isotropic TV, and SART-wavelet-FFT filter. A real sample was used to evaluate our algorithm, and all the algorithms were implemented in C++ language and matlab language on a PC (4.0-GB memory, 3.4-GHz CPU). Figure 1 shows a schematic diagram of an in-line phase contrast x-ray imaging setup. The projection data were collected from evenly distributed over pi-view in the BL13W1 of the Shanghai Synchrotron Radiation Facility (SSRF) in China. The images with a size of 1024 × 1024 pixels were reconstructed by SART, SART-isotropic TV, SART-wavelet-FFT filter, and SART-anisotropic TV. The relaxation parameter for SART was set to . The SART algorithm performed iterations K = 10.

3.1. The ring artifact correction method based on SART-isotropic TV

The TV-L2 model based on SART-isotropic TV in Eq. (22). We will illustrate some properties of the SART-isotropic TV based on numerical studies. The step size parameter β and the total number of gradient descents are chosen in order to obtain visually satisfactory results. We now discuss the impact of β for SART-isotropic TV. Particularly, we choose β from when we fix the total number of gradient descents.

As for SART-isotropic TV, the impact of β is illustrated in Fig. 3, in which we fix the total number of gradient descents . For the sake of better visualization, one region of interest (ROI) was magnified in Figs. 3(d)3(f), which shows that SART-isotropic TV is sensitive to the step size parameter β. We can see that large β leads to over-smoothed, and hence a moderate should be chosen for the best visualization.

Fig. 3. The performance of SART-isotropic TV with respect to β: (a) , (b) , (c) , and the corresponding magnified regions: (d) magnified ROI with respect to , (e) magnified ROI with respect to , and (f) magnified ROI with respect to .

We now discuss the impact of for SART-isotropic TV. Particularly, we choose from when we fix the step size parameter β = 0.02.

As for SART-isotropic TV, the impact of is illustrated in Fig. 4, in which we fix the β = 0.02. For the sake of better visualization, one region of interest (ROI) was magnified in Figs. 4(d)4(f), which shows that the ring artifact was suppressed as the increase of iteration index, and hence the should be chosen for the best visualization.

Fig. 4. The performance of SART-isotropic TV with respect to : (a) , (b) , (c) , and the corresponding magnified regions: (d) magnified ROI with respect to , (e) magnified ROI with respect to , (f) magnified ROI with respect to .
3.2. The ring artifact correction method based on SART-anisotropic TV

The TV-L2 model based on SART-anisotropic TV in Eq. (5). We will illustrate some properties of the SART-anisotropic TV based on numerical studies. SART-anisotropic TV is controlled via the regularization parameter λ. We choose different regularization parameter , , when we fix the maximum iteration number k = 110 of the solution of subproblems (11)–(15).

As for SART-anisotropic TV, the impact of λ is illustrated in Fig. 5, in which we fix the k = 110. For the sake of better visualization, one region of interest (ROI) was magnified in Figs. 5(d)5(f). Figure 5(e) shows that should be chosen for the balance between the ring artifacts reducing and the liver section details preserving.

Fig. 5. The performance of SART-anisotropic TV with respect to λ: (a) , (b) , (c) , and the corresponding magnified regions: (d) magnified ROI with respect to , (e) magnified ROI with respect to , and (f) magnified ROI with respect to .

We now discuss the impact of maximum iteration number k of the solution of subproblems (11)–(15) for SART-anisotropic TV. Particularly, we choose k from {110,140,170} when we fix the regularization parameter . As for SART-anisotropic TV, the impact of k is illustrated in Fig. 6.

Fig. 6. The performance of SART-anisotropic TV with respect to k: (a) k = 110, (b) k = 140, (c) k = 170, and the corresponding magnified regions: (d) magnified ROI with respect to k = 110, (e) magnified ROI with respect to k = 140, and (f) magnified ROI with respect to k = 170.

For the sake of better visualization, one region of interest (ROI) was magnified in Figs. 6(d)6(f), which shows that the SART-anisotropic TV succeeds in removing almost all visible rings to an appreciable level as the increase of iteration index. Figure 6(f) shows that the image is visually more satisfying, and hence the k = 170 should be chosen for the best visualization.

3.3. Comparison of ring artifact correction

First, we compare SART-anisotropic TV with two methods: SART and SART-isotropic TV. Figures 7(a)7(c) with a size of 1024 × 1024 pixels were reconstructed by different methods using 1024 projections. Figure 7(b) was reconstructed by SART-isotropic TV with β = 0.02 and . Figure 7(c) was reconstructed by SART-anisotropic TV with and k = 170.

Fig. 7. Reconstructed images using different methods for comparison (a) SART, (b) SART-isotropic TV, (c) SART-anisotropic TV, and the corresponding difference images for (d) the difference between SART and SART, (e) the difference between SART and isotropic TV, and (f) the difference between SART and anisotropic TV.

Figures 7(d)7(f) illustrate this effect by means of the difference images between images reconstructed by SART and methods: SART, SART-isotropic TV, and SART-anisotropic TV. Figure 7(e) shows that liver structures are visible in the difference image, which implies that SART-isotropic TV cannot preserve the liver section details. As can be seen from Fig. 7(f), no liver structures are visible in the difference image, which implies that SART-anisotropic TV can preserve spatial resolution.

Next, in order to provide further validation of the SART-anisotropic TV reconstruction accuracy, one region of interest (ROI) was magnified to allow a better detailed comparison. The magnified regions of the SART are illustrated in Fig. 8(d). The experiment result in Fig. 8(d) shows that the image reconstructed by SART is influenced by the ring artifacts. Figure 8(e) shows that ring artifacts correction conflicts with detail preservation for SART-isotropic TV. However, the results of reconstructions for the liver sample in Fig. 8(f) demonstrate that SART-anisotropic TV is effective and efficient for structure characteristics preservation on the ring artifacts correction process.

Fig. 8. Reconstructed images using different methods for comparison (a) SART, (b) SART-isotropic TV, (c) SART-anisotropic TV, and the corresponding magnified regions: (d) magnified ROI of SART, (e) magnified ROI of SART-isotropic TV, and (f) magnified ROI of SART-anisotropic TV.

Finally, in order to make a comparison study among the three algorithms: SART, SART-wavelet-FFT filter, and SART-anisotropic TV, two regions of interest(I and II)were magnified, which are shown in Fig. 9(a).

Fig. 9. (color online) Reconstructed images using different methods for comparison (a) SART, (b) SART-wavelet-FFT filter, (c) SART-anisotropic TV, and the corresponding magnified regions: (d) magnified ROI-I of SART, (e) magnified ROI-I of SART-wavelet-FFT filter, (f) magnified ROI-I of SART-anisotropic TV, (g) magnified ROI-II of SART, (h) magnified ROI-II of SART-wavelet-FFT filter, and (i) magnified ROI-II of SART-anisotropic TV.

As shown in Figs. 9(d) and 9(g), without sinogram preprocessing, the image reconstructed by SART is so strongly contaminated by ring artifacts, that in order to correct ring artifacts, the input sinogram is preprocessed using a wavelet-FFT filter. The wavelet-FFT filter algorithm,[7] in addition to the input sinogram, employs three parameters for the filtering process: the highest decomposition level L, the wavelet type, and the damping factor σ. The choice of these parameters permits the removal of most ring artifacts from the reconstructed slice, with minimal distortion at the image center. Higher order wavelets have a stronger effect on the center of the image, while decomposition levels smaller than 5 are instead not enough to remove all ring artifacts.[7] Taking this factor into account, the results have been obtained with the DB42 wavelet L =5 and σ = 2.4 (Fig. 9(b)). The experiment result in Fig. 9(h) shows that most ring artifacts are reduced from the reconstructed slice. However, Fig. 9(h) also shows that the image reconstructed by the SART-wavelet-FFT filter still contains a lot of noise. The experiment result in Fig. 9(e) shows that there are minimal distortions at the reconstructed image center(solid arrow). Compared with SART-wavelet-FFT filter, SART-anisotropic TV still gives a good indication of ring artifacts correction, which were shown in Figs. 9(f) and 9(i).

4. Conclusions

We proposed a new reconstruction method for IL-PC-CT. It is based on an anisotropic TV-L2 model to correct the ring artifacts and preserves fine structures during the image reconstruction process. The anisotropic TV regularization scheme ensures the proper balance between the extent of the ring reduction and the soft tissue detail. We applied SART and ADMM to solve the anisotropic TV-L2 model. The real experimental synchrotron data results show that the anisotropic TV-L2 model can improve the liver sample quality over the compared models in terms of preserving characteristics of the liver tissue and reducing the ring artifacts.

The ring artifacts can be of different types and of different intensities. On the one hand, completely damaged detector pixels cause strong ring artifacts in the tomographic image.[7] On the other hand, miscalibrated detector pixels, e.g., due to beam instabilities not completely taken into account by a normalization correction, give rise to wide and faint rings in the tomographic image.[34] In our study, the SART-anisotropic TV method is effective for the wide and faint ring artifacts. Further studies are currently under way to correct strong ring artifacts.

Reference
[1] Wang L Li X Wu M Zhang L Luo S Q 2013 BioMedical Engineering Online 12 1
[2] Boin M Haibel A 2006 Opt. Express 14 12071
[3] Herman G T Brouw W N 1980 Image Reconstruction from Projections — the Fundamentals of Computerized Tomography New York Academic Press
[4] Davis R Elliott J C 1997 Nucl. Instrum. Methods Phys. Res. 394 157
[5] Gorner W Hentschel M P Muller B R Riesemeier H Krumrey M Ulm G Diete W Klein U Frahm R 2001 Nucl. Instrum. Methods Phys. Res. 467 703
[6] Ketcham R A 2006 Proceedings of SPIE — The International Society for Optical Engineering 07 6318
[7] Münch B Trtik P Marone F Stampanoni M 2009 Opt. express 17 8567
[8] Sijbers J Postnov A 2004 Phys. Med. Biol. 49 147
[9] Prell D Kyriakou Y Kalender W A 2009 Phys. Med. Biol. 54 3881
[10] Zhang H M Wang L Y Yan B Li L Xi X Q Liu Y J 2013 Chin. Phys. 22 078701
[11] Fang S Guo H 2014 Chin. Phys. 23 057401
[12] Lin M B Zhao C X Yu X D Yu F S Yang T Y Yan K 2014 Acta Phys. Sin. 63 138701 in Chinese
[13] Gu Y F Yan B Li L Wei F Han Y Chen J 2014 Acta Phys. Sin. 63 018701 in Chinese
[14] Laroque S J Sidky E Y Pan X C 2008 J. Opt. Soc. Am. 25 1772
[15] Li S P Wang L Y Yan B Li L Liu Y J 2012 Chin. Phys. 21 108703
[16] Kostenko A Batenburg K J King A Offerman S E vanVliet L J 2013 Opt. Express 21 12185
[17] Chen Z Q Jin X Li L Wang G 2013 Phys. Med. Biol. 58 2119
[18] Andersen H Kak A C 1984 Ultrason. Imaging 6 81
[19] Wang D Zeng X J Keane J A 2006 Danyluk Et 15 92
[20] Rudin Leonid I Osher S Fatemi E 1992 Physica D: Nonlinear Phenomena 60 259
[21] Esedoglu S Osher S J 2004 Commun. Pure Appl. Math. 57 1609
[22] Lou Y Zeng T Osher S Xin J 2015 SIAM J. Imag. Sci. 8 1798
[23] Torres-Forné A, Marquina A, Font J A and Ibánez J M 2016. arXiv: 1602.06833[astro-ph.IM]
[24] Kang R R Qu G R Wang B 2016 Circuits Systems Signal Process 35 3380
[25] Vandeghinste B Goossens B Van H R Vanhove C Pizurica A Vandenberghe S Staelens S 2013 IEEE Trans. Nucl. Sci. 60 3305
[26] Choksi R van Gennip Y Oberman A 2011 Inverse Problems Imaging 5 591
[27] Boyd S Parikh N Chu E Peleato B Eckstein J 2011 Foundations and Trends in Machine Learning 3 1
[28] Jian J Yang H Zhao X Xuan R Zhang Y Li D Hu C 2013 J. Synch. Radiat. 23 600
[29] Parikh N Boyd S 2013 Foundations and Trends in Optimization 1 123
[30] Gao H 2016 Phys. Med. Biol. 61 7187
[31] Liu L K Chan S Nguyen T 2015 IEEE Trans. Image Process. 24 3905
[32] Ng M K Chan R H Tang W C 1999 SIAM J. Sci. Comput. 21 851
[33] Li C 2009 An efficient algorithm for total variation regularization with applications to the single pixel camera and compressive sensing Ph. D. dissertation Rice University
[34] Raven C 1998 Rev. Sci. Instrum. 69 2978